#this programs reads in the composite dose dicom RTDose file and replaces the pixel array 
#dose values with BED values calculated in Mathematica
import os #imports the operating system module, needed for file and directory functions
import numpy as np #imports the NumPy scientific computing library and makes np shorthand for numpy
os.chdir("/Users/Documents/pydicom-0.9.7") #sets the directory path
import dicom #imports the pydicom package
composite_dose_file=dicom.read_file("composite dose.dcm") #reads in the composite dose dicom file
#the line below reads in the VoxelBedsForDicom text file produced in Mathematica
beds_text_file_import=np.fromfile('VoxelBedsForDicom.txt', dtype=int,count=-1,sep=" ")
#the line below reshapes the one-dimensional text file array into an 11x10x9 array
beds_array=np.array(beds_text_file_import).reshape((11, 10,9)) 
max_bed=np.amax(beds_array) #finds the maximum array value
#the line below divides the maximum allowed array value (65535) by max_bed
#(decimal point needed for floating point division in Python 2)
array_scalingfactor=65535./max_bed
#the line below scales the array values to prevent overflow of values over 65535
adjusted_beds_array=beds_array*array_scalingfactor 
#the for loop below iterates through the array values rounding and converting them to integers
#and replaces the values in composite_dose.pixel_array with these new values
for i in range(0,11):
    for j in range(0,10):
        for k in range(0,9):
            composite_dose_file.pixel_array[i,j,k]=int(round(adjusted_beds_array[i,j,k]))
#the line below creates a new Dose Grid Scaling factor adjusting the original Dose Grid Scaling 
#factor with the array_scalingfactor
scalingfactor=float(composite_dose_file[0x3004,0x0e].value)/array_scalingfactor
#the line below rounds scalingfactor to 11 decimal places, converts it to a string
#and sets the Dose Grid Scaling factor to this value
composite_dose_file[0x3004,0x0e].value=str(round(scalingfactor,11))
#the line below overwrites the PixelData attribute of the composite_dose_file changing the original
#dose values to BED values
composite_dose_file.PixelData=composite_dose_file.pixel_array.tostring()
#the line below saves the modified "composite dose" dicom file as "composite beds"
#and thus the original RTDose file has been updated with BED values rather than dose values
composite_dose_file.save_as("composite beds.dcm")